######################################################################################
######################################################################################
#Replication code for models presented in:

#Buntaine, Mark T. and Bradley C. Parks. 
#2013. 
#When Do Environmentally Focused Assistance Projects Achieve their Objectives? Evidence from World Bank Post-Project Evaluations. 
#Global Environmental Politics 13(2): 67-90.

#Mark Buntaine contact (as of March 2013): mbuntaine@wm.edu / 919.627.7066
#Bradley Parks contact (as of March 2013): bcpark@wm.edu / 

#Current for R Version 2.15.2 (version "Trick or Treat")
######################################################################################
######################################################################################

##########################################################################################
###Setting Up
##########################################################################################

library(arm) #Version 1.5-08
#Load "GEP_2013_replication_data.csv" into R as object named "data"

data$global.percent <- data$global.prop/100
data$tri.outcome <- as.factor(data$tri.outcome)

##########################################################################################
###TABLE 2
##########################################################################################

data.mod1 <- na.omit(subset(data, select=c("tri.outcome","WBI.govteff","global.prop","IEG.qs.tri","env.damage","env_pct","concessional","FH.civlib","ENGOs.pmil","amount.mil","IEG.EvalType")))
mod1 <- polr(tri.outcome ~ WBI.govteff + IEG.qs.tri + env_pct + global.prop + env.damage + concessional + FH.civlib + ENGOs.pmil + amount.mil + IEG.EvalType, data.mod1, Hess=T)
summary(mod1)
#Here is example code to get p-value: 2*pt(-1.92,df=mod1$df.residual)

#Example of how to calculate null deviance:
mod1.null <- polr(tri.outcome ~ 1, data.mod1, Hess=T)
summary(mod1.null)

mod2<-step(mod1,direction="both")
summary(mod2)

#Code to examine correlate between indepedent variables
cor <- data.mod1[,2:10]
cor(cor,use="pairwise.complete.obs")
#Removing: concessional, FH.civlib for mod3 as described in the paper

data.mod3 <- na.omit(subset(data, select=c("tri.outcome","WBI.govteff","global.prop","IEG.qs.tri","env.damage","env_pct","ENGOs.pmil","amount.mil","IEG.EvalType")))
mod3 <- polr(tri.outcome ~ WBI.govteff + IEG.qs.tri + env_pct + global.prop + env.damage + ENGOs.pmil + amount.mil + IEG.EvalType, data.mod3, Hess=T)
summary(mod3)

data.mod4 <- na.omit(subset(data, select=c("tri.outcome","WBI.govteff","global.prop","IEG.qs.tri","env_pct","amount.mil","IEG.EvalType")))
mod4 <- polr(tri.outcome ~ WBI.govteff + IEG.qs.tri + env_pct + global.prop + amount.mil + IEG.EvalType, data.mod4, Hess=T)
summary(mod4)

##########################################################################################
###TABLE 3
##########################################################################################

data.mod5 <- na.omit(subset(data, select=c("IEG.outcome.bin","IEG.outcome.bin2","WBI.govteff","global.prop","IEG.qs.tri","env.damage","env_pct","concessional","FH.civlib","ENGOs.pmil","amount.mil","IEG.EvalType")))

mod5a <- glm(IEG.outcome.bin ~ WBI.govteff + IEG.qs.tri + env_pct + global.prop + env.damage + concessional + FH.civlib + ENGOs.pmil + amount.mil + IEG.EvalType, family=binomial, data=data.mod5)
summary(mod5a)

mod5b <- glm(IEG.outcome.bin2 ~ WBI.govteff + IEG.qs.tri + env_pct + global.prop + env.damage + concessional + FH.civlib + ENGOs.pmil + amount.mil + IEG.EvalType, family=binomial, data=data.mod5)
summary(mod5b)

#Step-AIC on the models
mod6a <- step(mod5a, direction="both")
summary(mod6a)
mod6b <- step(mod5b, direction="both")
summary(mod6b)

data.mod7 <- na.omit(subset(data, select=c("IEG.outcome.bin","IEG.outcome.bin2","WBI.govteff","global.prop","IEG.qs.tri","env_pct","concessional","amount.mil","IEG.EvalType")))

mod7a <- glm(IEG.outcome.bin ~ WBI.govteff + IEG.qs.tri + env_pct + global.prop + concessional + amount.mil + IEG.EvalType, family=binomial, data=data.mod7)
summary(mod7a)

mod7b <- glm(IEG.outcome.bin2 ~ WBI.govteff + IEG.qs.tri + env_pct + global.prop + concessional + amount.mil + IEG.EvalType, family=binomial, data=data.mod7)
summary(mod7b)